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Experimental measurements and numerical simulations are made on fluid flow through cross-slot 
devices with a range of aspect (depth:width) ratios, 0.4 < Q < 3.87. For low Reynolds numbers Re, 
the flow is symmetric and a sharp boundary exists between fluid streams entering the cross-slot from 
opposite directions. Above an a-dependent critical value 20 < Rec(a) < 100, the flow undergoes a 
symmetry-breaking bifurcation (though remains steady and laminar) and a spiral vortex structure 
develops about the central axis of the outflow channel. An order parameter characterizing the 
instability grows according to a sixth-order Landau potential, and shows a progression from second 
order to first order transitions as a increases. A tricritical point occurs for a « 0.55. The spiral 
vortex acts as a mixing region in the flow field and this phenomenon can be used to drive enhanced 
mixing in microfluidic devices. 

PACS numbers: 05.70.Fh, 47.20.Ky, 47.61.Ne, 64.60.Kw 


The ability of fluids to mix is greatly enhanced by tur¬ 
bulence, which occurs at large values of the Reynolds 
number Re = ULjv, where U and L are characteristic 
velocity and length scales respectively and v is the kine¬ 
matic viscosity of the fluid. Small length scales tend to 
suppress Re, making it difficult to develop turbulent mix¬ 
ing in microfluidic devices. Achieving efficient mixing via 
diffusion dominated processes presents a major techno¬ 
logical challenge to the expanding field of lab-on-a-chip 
development. 


The cross-slot device is emerging as a promising plat¬ 
form for enhancing mixing of fluids in micro-scale geome¬ 
tries 0,0. The planar elongational flow field generated 
by the cross-slot geometry (Fig. 1(a)) has found appli¬ 
cations in many research areas including for studies of 
macromolecular dynamics extensional rheometry 

and elastic instabilities of viscoelastic fluids , hydro- 
dynamic trapping 0, and imposing controlled defor¬ 
mations to complex biological structures [lll - [l^ . It has 
been known since the early 1990’s that such intersecting 
flows are prone to instability beyond a modest critical 
Reynolds number Rec ~ 0(10-100). In very deep cross¬ 
slot flow channels the observed instability takes the form 
of a stack of three-dimensional vortical structures that 
appear in the central cross-over region 1^, [16| , while in 
the related 4-roll mill apparatus a flow field incorporat¬ 
ing a single helical region has been observed [3. In this 
work, with cross-slots of modest aspect ratios typical of 
microfluidic devices (depth to width ratios between 0.4 
and 3.87), we have found that the instability occurs above 
an aspect ratio-dependent value of Rec and results in a 
simple straight vortex that extends downstream along the 
outlet channels, see Fig. [ijb). This flow instability has 
been shown to promote mechanical scission of polymer 
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FIG. 1; (Color online) (a) Schematic diagram of a cross-slot 
device. Flow enters along y and exits along x. Confocal mi¬ 
croscopy is performed in z-planes, which are scanned through 
the full depth of the device and used to reconstruct an im¬ 
age in the x = 0 plane (shaded (green) region), (b) Three- 
dimensional (3D) rendering of a vortex structure observed for 
the flow of water at Re = 75.8 in a cross-slot with a = 1. The 
image is generated from 2 -plane images spaced at Sz = 5pm 
and has been cropped around the central vortex. The vol¬ 
ume shown corresponds to the fluorescently-dyed fluid stream. 
Also see Movie Ml for an animated version of Fig. 1(b) p^ . 




chains |18l| , and has more recently been demonstrated to 
enhance mixing between the two incoming fluid streams 
in micro-sized devices 0 , 0 . Numerical simulations have 
shown that the mixing quality is higher and Rec is lower 
in the cross-slots than in the more well-known T-shaped 
mixing device with equivalent channel dimensions 0. 

In this Letter, we report the results of detailed exper¬ 
imental and numerical studies of the spiral vortex flow 
instability in cross-slots with a range of aspect ratios and 
over a wide range of Re. In contrast to previous stud¬ 
ies 0E [l6l |. we identify appropriate order parameters 
that characterize the instability as a function of Re in 
each case. In particular, we present a systematic analysis 
in terms of bifurcation theory analogous to the Landau 
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theory of phase transitions. The observed phenomena 
are well described by a Landau-type sixth-order polyno¬ 
mial potential 23 - 22 |, with parameters that show the 
transition develops from second order to first order as a 
increases, passing through a tricritical point for a ~ 0.55. 
These data can be fully described by scaling theory, and 
the universal scaling function and behavior are measured 
near the tricritical point. Improved understanding and 
characterization of stability conditions for flows through 
intersecting geometries is vital for the optimization of 
many laboratory microfluidic experiments and also prac¬ 
tical lab-on-a-chip designs, including for the specific goal 
of enhancing the mixing of fluids in channels with small 
dimensions operating at low Re. 

The cross-slot device (Fig.[T][a)) consists of two bisect¬ 
ing rectangular channels of width w and depth d, with the 
aspect ratio defined by a = d/w. In the experiments four 
devices are utilized, all with d = 1.2 mm but with w var¬ 
ied in order to provide a = 0.49,1.00,1.85,andS.87. An 
inlet length of at least 12.5w ensures a fully-developed 
flow before the fluid reaches the central region of each 
device. Newtonian fluid (water) is pumped into two op¬ 
posing channels (along the y-direction) and exits through 
the two opposing outlet channels (along the x direction). 
One of the incoming fluid streams is fluorescently-dyed 
with rhodamine B (concentration = 10 /iM, Sigma) and 
a laser-scanning confocal microscope (Zeiss LSM 780) is 
employed to examine the interface between fluid streams 
where they meet in the central cross-over region. Imag¬ 
ing in closely-spaced planes through the depth of the 
flow cell (—(i/2 < z < d/2) allows accurate recon¬ 
struction of an image in the a; = 0 plane, see (green) 
shaded area in Fig. [TJa). Experiments are performed 
at 24°C over a range of Re = pwU/fx, with the fluid 
density p = 997.1 kgm“3 and the dynamic viscosity 
pL = 9.1 X 10“^ Pas. The average flow velocity U within 
each channel of the cross-slot is controlled using a pre¬ 
cision dual syringe pump (Harvard PHD Ultra). While 
for low Re the interface between fluid streams is sharp 
and vertical over the y = 0 plane, beyond a fairly mod¬ 
erate critical value RCc the flow bifurcates and breaks 
symmetry (but remains steady and laminar) and intri¬ 
cate spiral vortex structures develop, see Fig. [TJb) and 
Movie Ml 0. 

The numerical method solves the equations of motion 
pDu/Dt = — Vp-f/iV^u and mass conservation V- u = 0, 
for laminar flow of a Newtonian incompressible fluid by 
using a fully-implicit, second-order finite volume method 
0, liS Hi- For the inertial flows considered here, rel¬ 
evant features of the method are the treatment of the 
partial time derivative du/dt using the three time-level 
pressure correction algorithm and the treatment of 
the convective term u- Vu, which is discretized using the 
3rd order high-resolution CUBISTA scheme [i^. In this 
way we maintain 2nd-order accuracy in space and time. 

Non-uniform orthogonal meshes are deployed on the 


3D cross-slot geometry of Fig. 1(a), with the central cube 
(for a = 1) having 25^ uniform control volumes (CV’s) 
on the base mesh, and 35^ or 71^ CV’s on more refined 
meshes. The CV size expands slowly towards the two 
inlets, along the channel aligned with the y-axis, and to¬ 
wards the two outlets, along the x-axis (see Fig. SI 0 ). 
Theoretical [23 2D fully-developed velocity {vth{x,z)) 
profiles are applied at the inlets, while at the outlets 
zero axial gradients are assumed for velocity components 
{d/dx = 0), pressure is linearly extrapolated and the to¬ 
tal flow rate is forced to satisfy overall mass conservation 
{Qout = ‘2Qin, where Qin = Uwd for each inlet section). 
Note that we do not force the flow rate entering via each 
inlet to divide equally between each outlet, but only that 
the flow rate through each outlet is equal, as in the exper¬ 
iments. The numerical simulations explore devices with 
aspect ratios set equal to the four experimental devices 
as well as additional values of aspect ratio near to the 
tricritical value. 

Great care is exercised to guarantee that the numerical 
results are well converged and do not significantly with 
the computational mesh. In Fig. SI [l9|, we illustrate 
the various meshes employed and compare the results ob¬ 
tained for a cross-slot with a = 1, showing good accuracy. 

In Fig. [5Ka-d) we present reconstructed confocal im¬ 
ages of the a; = 0 plane (spanning —w/2 < y < w/2, 
—d/2 < z < d/2) for the cross-slot device with a = 1, 
which show the evolution of the vortical structure as Re 
increases. Qualitatively similar behavior is observed in 
all four geometries (see Fig. S2 and Movies M2-M5 0): 
at low Re the interface between fluid streams is sharp 
and vertical as expected, however as Re increases above 
a critical value the spiral vortex forms abruptly about the 
x-axis. With further increases in Re, the spatial extent 
of the spiral expands, as does the number of turns of the 
arms. The value of the critical Reynolds number llMc,exp 
for the appearance of the spiral in the experiments de¬ 
pends on the aspect ratio, but is broadly consistent with 
values reported previously in related experiments iEl- 
17| (Rec,ea:p ~ 100, 40, 24 and 26 for a = 0.49, 1.00, 


1.85 and 3.87, respectively). In each of the four devices 
the central vortex forms with a favored orientation about 
the x-axis (anticlockwise for a = 0.49 and a = 1.00 and 
clockwise for a = 1.85 and a = 3.87, see Fig. S2). The 
orientation is presumably biased by some minor geomet¬ 
rical imperfections in the devices and causes the bifurca¬ 
tion to follow a favored branch. In only a few instances 
in the cross-slots with a = 1.85 and a = 3.87 at higher 
Re did we observe the vortex to form in the unfavored 
sense (see Fig. S3 0 for examples). The vortex struc¬ 
tures observed in the a; = 0 plane can be well-fitted by 
elliptical logarithmic spirals (see Fig. S4 0). 

Numerical simulations result in remarkably good 
agreement with the experiment, as shown by the stream¬ 
line plots in Fig. [Uc-h) that can be directly compared 
with the experimental results in Fig. [2Ka-d). Note that 
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FIG. 2: (Color online) Confocal imaging (a-d) and numerically generated streamlines (e-h) depicting the evolution of flow 
structures in the x = 0 plane for Newtonian fluid flow in the cross-slot device with a = 1 under the following conditions: (a, 
e) Re = 15.2, (b, f) Re = 42.8, (c, g) Re = 60.6, and (d, h) Re = 91.0. In (a-d) Fluorescently-dyed fluid enters from the right 
(positive y) and undyed fluid enters from the left (negative y); outflow is along x (i.e. normal to the page). Scale bar in (a) 
represents 200 ^m. In (e-h) the streamlines are colored in order to resemble the experiment. The numerical result shown in 
(f) is one of two possible solutions at this Re; a symmetric solution can also be obtained, as shown in Fig. S5 [l^ . (i) Mixing 
parameter in all four experimental cross-slot devices, evaluated over the region spanned by —w/2 < y < w/2, —d/2 < z < d/2 
in the a: = 0 plane according to Eq. (1). (j) Numerical order parameter as a function of the control parameter, e, fitted with a 
Landau sixth order polynomial potential Eqs. (3, 4); broken lines represent unstable branches. 


in the numerical simulations at a = I.OO there is hys¬ 
teresis in the transition and both symmetric and asym¬ 
metric steady solutions can be obtained for Re = 42.8 
(Fig. [2If)), depending on whether Re is quasistatically 
increased or decreased from below or above the transi¬ 
tion, respectively. Both possible solutions are presented 
in Fig. S5 [l^. Additional numerical plots of velocity 
vector fields at various Re spanning the transition for the 
case of a = 1.00 are provided in Fig. S6 [l^. Our experi¬ 
mental method is not amenable to a genuinely quasistatic 
variation of Re since the flow is necessarily interrupted 
between image acquisition at each Re in order to refill 
the syringes with fluid. Therefore our experiment can 
not resolve the hysteresis in the transition. 


We treat our experimental data following the approach 
suggested by Stroock et al. [2^ in order to assess the de¬ 
gree of mixing between the incoming fluid streams, M^^p- 
This parameter is computed using the standard deviation 
of the pixel intensity evaluated over images such as those 
shown in Fig. [21 a-d) and Fig. S2 [l^ : 

Me.p = l-2((/-(/))2)l/2, (1) 

where I is the grayscale pixel value (normalized between 
0 and 1, corresponding to the minimum and maximum 
grayscale intensities in the images) and () indicates an 
average over all pixels in the field of view. For com¬ 
pletely segregated fluid streams, there is a binary distri¬ 
bution of pixel intensities, the standard deviation = 0.5 
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and hence M^xp = 0; this is essentially the case illus¬ 
trated in Fig. IHa). Complete mixing between the fluid 
streams would mean a uniform pixel intensity of / = 0.5 
over the entire field of view, therefore a standard devia¬ 
tion of zero and hence M^xp = 1- In Fig.[2fi), we plot the 
mixing parameter, M^xp, as a function of the Reynolds 
number for all four experimental aspect ratios. Fig. mi) 
shows a general reduction in Ksc^exp and an increasingly 
large and rapid increase in as a increases. We also 
note that the transition is imperfect; there is a small but 
distinct curvature in M^xp near the onset, which is con¬ 
sistent with the presence of their being a favored branch 
for the transition (clockwise or anticlockwise spirals). We 
treat the favored branch as the positive branch regardless 
of the spiral orientation. 

In the case of the numerical simulations, we find the 
growth of the instability is best described by the following 
parameter: 

M,.„„ = 1;=^, (2) 

where Vmax\x=y=o maximum value of the y- 

component of the velocity measured along the z-axis at 
each Re (for a symmetric flow Mnum = 0). A pictorial 
description of Mnum is given in Fig. S6 [1^. Fig. [2Kj) 
shows a plot of Mnum as a function of the control param¬ 
eter £ = (Re—Reg)/Reg for various aspect ratio cross-slot 
devices. The data in Fig. [2Kj) is fitted using a Landau- 
type model with ’’free energy” F, given by: 

F = -htl) - -f 1?/’^ -b (3) 

where the order parameter tjj = Mnum-, in the numerical 
case. Even terms are included in Eq. (3) since for a per¬ 
fect system F should be independent of the sign of ip (i.e. 
the handedness of the spiral). The lowest order asym¬ 
metric “field” term [hip) can account for imperfections 
that bias the handedness. In the numerical simulations 
the vortex can form in the clockwise or anticlockwise ori¬ 
entation with similar probablity, hence h = 0 (for all 
aspect ratios). For given values of the parameters h, g, 
and fc, the value of ip{e) corresponds to the extrema of 
F, dF/dip = 0, giving: 

£ = Re/ROg — 1 = kip^ + gip'^ — hip~^. (4) 

The ratio of the coefficients g/k > 0 ior a Pp, 0.55 (cor¬ 
responding to a forward bifurcation or a second-order 
transition) and decreases monotonically with increasing 
a, turning negative for a > 0.55 at which point the bi¬ 
furcation becomes backwards, or first order and an in¬ 
creasing large hysteresis loop grows as a increases. The 
numerical data obtained for a = 0.55 are well-fitted by 
Eq. (4) with g = 0, (cyan) stars in Fig.[2Kj), and thus cor¬ 
responds to a tricritical transition. Fitting of the numer¬ 
ical data with Eq. (4) provides the values of RCg shown 




r = ek/g^ 

FIG. 3: (Color online) (a) Experimental mixing parameter 
fitted with the sextic Landau potential indicated by the nu¬ 
merical results, (b) Collapse of experimental data (solid sym¬ 
bols) and numerical data (hollow symbols) onto the theoreti¬ 
cal curve given by Eq. (5) (solid lines). 

in Fig. S7 [l^. For a > 0.55, due to the hysteresis loop 
we can find two values of the critical Reynolds number: 
Reg denotes the value found for quasistatic increases in 
Re from below, while Re* denotes the value found for 
quasistatic decreases in Re from above the onset. The 
numerically predicted hysteresis loops for a ^ 0.55 can 
not be observed experimentally, however we find that our 
experimental critical Reynolds numbers (Rec^g:jp) agree 
very well (within approx. ±5%) with the numerical val¬ 
ues of Re* (Fig. S7). 

In Fig. Ela) we fit our experimental mixing parame¬ 
ter with Eq. (4), setting the order parameter ip = Mf^xp 
and fixing Reg to the value determined numerically for 
each aspect ratio device. We introduce a small coeffi¬ 
cient 0.002 < h < 0.003 into Eq. (4) that accounts for 
geometrical imperfections in the experimental devices. In 
general we find that our experimental order parameter is 
well described by the same sextic Landau potential that 
describes the numerical data. 

The theory for tricritical points predicts a univer¬ 
sal scaling form for the canonical order parameter (p = 
ip\/k/\g\ in terms of the control parameter r = ekjg^, 
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showing data collapse in the form: 

cj/ = {l- sgn( 5 )Vl + 4r)/2. (5) 


In Fig.[31Jb), we have plotted both the experimental and 
numerical results in scaled form and show the compar¬ 
ison with the data collapse prediction. The agreement 
is excellent and confirms our identification of the mixing 
transition as a tricritical one. 

In summary, we have demonstrated that the spiral in¬ 
stability observed for Newtonian flow in the cross-slot de¬ 
vice driven far from equilibrium above a critical Reynolds 
number, can be well-described by a Landau model anal¬ 
ogous to that used near equilibrium tricritical points. 
We note a certain similarity of our results with those 
presented for flow in low-aspect ratio Taylor-Couette de¬ 
vices above a critical angular velocity (e.g. 21, 22, S^), 
however the instability in the cross-slot geometry can be 
harnessed to promote mixing of Newtonian fluids at the 
quite modest Reynolds numbers accessible in microflu¬ 
idic devices. Such improved mixing could benefit the 
efficiency of lab-on-a-chip applications such as chemical 
synthesis, antibody antigen binding reactions and bioas¬ 
say. It is likely that bifurcation phenomena reported 
in related flow geometries (such as in the T-shaped mi¬ 
cromixer, for example [2,1^) could also be characterized 
by the tricritical point formalism presented here. 
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